home *** CD-ROM | disk | FTP | other *** search
/ Power Programmierung / Power-Programmierung (Tewi)(1994).iso / magazine / progjour / 1988 / 06 / fcode / ssu.for < prev    next >
Text File  |  1987-09-03  |  13KB  |  432 lines
  1. C./ ADD NAME=SSU
  2. C     SUBROUTINE TO TEST ACCURACY OF ROOTS OF CUBIC EQUATION,
  3. C
  4. C     SUBROUTINE ROOTST ( K, KAPPA, H )
  5. C     IMPLICIT REAL*8 (A-H,O-Z)
  6. C     REAL*8 K, K2, K3, KAPPA
  7. C
  8. C     K2 = K*K
  9. C     K3 = K2*K
  10. C     T1 = DABS(2.*K3)
  11. C     T2 = H*K2
  12. C     DT = DABS(T1-DABS(KAPPA-T2))
  13. C     T = T1
  14. C     IF (DABS(KAPPA).GT.T) T = DABS(KAPPA)
  15. C     IF (DABS(T2).GT.T) T = DABS(T2)
  16. C     DT = DT/T
  17. C     IF (DT.LT..0001) RETURN
  18. C     PRINT 1000, KAPPA, H, K, DT
  19. C1000 FORMAT (' ERROR IN ROOT, KAPPA = ',E12.4,' H = ',E12.4,'
  20. C    * K = ',E12.4,' DT = ',E12.4 )
  21. C     END
  22. C
  23. C     SUBROUTINE TO CALCULATE THE COMPLETE ELLIPTIC INTEGRAL
  24. C     OF THE THIRD KIND,
  25. C
  26.       SUBROUTINE CEL3( CE, KP, M, FP )
  27.       IMPLICIT REAL*8 (A-H,O-Z)
  28.       REAL*8 KC, M, M0, KP
  29.       KC = KP
  30.       P = FP
  31.       IF (KC*P.EQ.0.D0) GO TO 10
  32.       KC = DABS(KC)
  33.       E = KC
  34.       M0 = 1.D0
  35.       IF (P.LE.0.D0) GO TO 2
  36.       C = 1.D0
  37.       P = DSQRT(P)
  38.       D = 1.D0/P
  39.       GO TO 3
  40.  2    G = 1.D0 -P
  41.       F = KC*KC - P
  42.       P = DSQRT(F/G)
  43.       D = -M/(G*P)
  44.       C = 0.D0
  45.  3    F = C
  46.       C = C + D/P
  47.       G = E/P
  48.       D = 2.D0 * ( F*G + D )
  49.       P = G + P
  50.       G = M0
  51.       M0 = KC + M0
  52.       IF ( DABS(G-KC).LE.1.D-4*G ) GO TO 4
  53.       KC = 2.D0 * DSQRT(E)
  54.       E = KC * M0
  55.       GO TO 3
  56.  4    CE = 1.5707963D0 * (C*M0+D)/(M0*(M0+P))
  57.       RETURN
  58.  10   CE = (1.D18)**2
  59.       RETURN
  60.       END
  61. C
  62.       SUBROUTINE SS( E, DE, S, DELTA, EC, BET )
  63. C     SIGNALS
  64. C     DLE = .0 AND DE NOT= 0, NORMAL,
  65. C     DLE > .0, COMPUTES COEFFICIENT OF LN(E-EL) NEAR E = EL,
  66. C     DLE < .0, COMPUTES COEFFICIENT OF (E-EL)**1/3 FOR E NEAR
  67. C     EL AND EL = EMID,
  68. C     DE = .0, GIVES END CORRECTIONS, S FOR E < OR E > EMAX,
  69. C
  70.       IMPLICIT REAL*8 (A-H,O-Z)
  71.       real*8 darcos
  72.       DIMENSION S(3), X(3)
  73.       REAL*8 KAPPA, KAPRM, K(3), KL, NUM(3), K2, KC, KK, KC2
  74.       COMMON KAPPA, KAPRM, H, Q, RQ, RP, H2, R1
  75.       COMMON/A/NA, NS, KL, ISIG
  76.       COMMON/C/EPS, DLE
  77.       COMMON/D/S12, S11 /PI/PI
  78.       PI = DATAN(1.D0)*4.
  79.       DO 1 N = 1, 3
  80.  1    S(N) = 0.D0
  81.       DDE = .1*DE
  82.       E2 = EC - DELTA/3.D0
  83.       E3 = EC + DELTA/3.D0
  84.       EMIN = .0
  85.       IF (E2.LT.0.D0) EMIN = E2
  86.       EMAX = .0
  87.       IF (E3.GT.0.D0) EMAX = E3
  88.       EMID = E2
  89.       IF (EMID.EQ.EMIN) EMID = .0
  90.       IF (EMID.EQ.EMAX) EMID = E3
  91.       IF (DABS(E).LT.DDE) E = .0
  92.       IF (DABS(E-E2).LT.DDE) E = E2
  93.       IF (DABS(E-E3).LT.DDE) E = E3
  94.       B2 = BET**.66666666D0
  95.       B4 = B2**2
  96.       RP = B4 + 2.D0/B2
  97.       RQ = 18.D0/RP
  98.       IF (DLE.GT.0.) E = E-EPS
  99.       H = RP*E-2.D0*EC/B2
  100.       H2 = 6.D0*E - 4.D0*EC
  101.       Q = B4*(2.D0*H**2/3.D0)-BET**2*E**2-2.*((E-EC)**2+DELTA**2/9.D0)
  102.       Q = 3.D0*Q/(RP*B4)
  103.       R1 = 2.D0*H*(2.*BET**2+1.D0)/(RP*B2)
  104.       KAPPA = E*(E-E2)*(E-E3)
  105.       KAPRM = E*(E-E2) +E*(E-E3) +(E-E2)*(E-E3)
  106.       XP = KAPRM/RP
  107.       Y = ( Q + RQ*XP )
  108.       YP = ( R1 + 18.D0*H2/RP**2)
  109.       P23 = 2.D0*3.1415927D0/3.D0
  110.       IF (DLE.GT.0.D0) GO TO 50
  111. C     IF (DLE.EQ.0.D0) GO TO 7
  112. C     EL = EMID
  113. C     E = EL
  114. C     XA = ((EMID-EMIN)*(EMAX-EMID)*DE/2.D0)**.66666666D0
  115. C     A = DSQRT(3.D0)*XA
  116. C     RA = DSQRT(A)
  117. C     AK = DSQRT((1.D0-DSQRT(3.D0)/2.)/2.)
  118. C     AP = A+XA-XP
  119. C     CALL CEL1(KK, AK, IER)
  120. C     FJ1 = -KK/RA
  121. C     FJ2 = FJ1/(XP*RP)
  122. C     FJ3 = FJ2/(XP*RP)
  123. C     EK = .0
  124. C     GO TO 30
  125. C
  126. C     NORMAL CASE
  127.  7    G = H/3.D0
  128.       IF (DE.EQ.0.) GO TO 40
  129.       IF (DABS(H**3).LT.DABS(KAPPA*1.E-10)) GO TO 4
  130. C     THREE ROOT REGION
  131.       C = 2.D0* KAPPA*(3.D0/H)**3 -1.D0
  132.       IF (C*C.GE.1.) GO TO 10
  133.       PHI = DARCOS(C)/3.D0
  134.       K(1) = DABS( G*(DCOS(PHI) -.5D0 ))
  135. C
  136.       K(2) = DABS( G*(DCOS(PHI-P23) -.5D0 ))
  137. C
  138.       K(3) = DABS( G*(DCOS(PHI+P23) -.5D0 ))
  139. C
  140.       IF (ISIG.GT.1) PRINT 1000, (K(L), L = 1, 3 )
  141.  1000 FORMAT (' K S = ', 3E12.4 )
  142.       DO 3 L = 1,3
  143.  3    X(L) = K(L)**2
  144.       IF(DABS(E).LT.DDE.OR.DABS(E-E2).LT.DDE.OR.DABS(E-E3).LT.DDE)
  145.      *  GO TO 18
  146.       K2 = (X(2)-X(1))/(X(3)-X(1))
  147.       AK = DSQRT(K2)
  148.       ALPHA2 = (X(2)-X(1))/(XP-X(1))
  149.       KC2 = (X(3)-X(2))/(X(3)-X(1))
  150.       KC = DSQRT(KC2)
  151.       FP = (XP-X(2))/(XP-X(1))
  152.       CALL CEL1( KK, AK, IER )
  153.       CALL CEL2( EK, AK, 1.D0, KC2, IER )
  154.       CALL CEL3( PIK, KC, K2, FP )
  155.       W = ( KAPPA - H*XP )**2 - 4.*XP**3
  156.       RTW = DSQRT(W)
  157.       WP = -2.*H2*(6.D0*XP**2-XP*H**2+H*KAPPA)/RP
  158.       RAC = DSQRT(X(3)-X(1))
  159.       X1P = X(1) -XP
  160.       X2P = X(2) -XP
  161.       X3P = X(3) -XP
  162.       FJ1 = -2.*KK/RAC
  163.       FJ2 = 2.*PIK/(RP*RAC*X1P) - PI/(RP*RTW)
  164.       FJ3 = -2.*RAC*EK/(W*RP**2)+.5D0*KK/(RAC*X1P*X2P*RP**2)
  165.      * -(6.D0*XP**2-XP*H**2+H*KAPPA)*PIK/(RAC*W*X1P*RP**2)
  166.       FJ3 = 2.*FJ3+PI*(6.D0*XP**2-XP*H**2+H*KAPPA)/(W*RTW*RP**2)
  167.       S(1) = YP*FJ2-H2*Y*FJ3
  168.       S(2) = 18.D0*FJ1/RP**2-Y*FJ2
  169.       S(3) = 2.*(X3P*KK/RAC-RAC*EK)*RP**2/3.D0
  170.       GO TO 20
  171. C
  172. C     ONE ROOT REGION
  173.  4    K(3) = (DABS(KAPPA/2.))**.33333333D0
  174. C
  175.       GO TO 11
  176.  10   TEMP = ( DABS(C) + DSQRT(C*C-1.D0))**.33333333D0
  177.       TEMP = DSIGN(TEMP,C)
  178.       K(3) = DABS( .5D0*( TEMP -1.D0 +1.D0/TEMP )*G)
  179. C
  180.  11   K(1) = -1.
  181.       IF (ISIG.GT.1) PRINT 2000, K(3)
  182.  2000 FORMAT (' K(3) = ',E12.4 )
  183.  18   XA = K(3)**2
  184.       A2 = (6.D0*XA**2-XA*H**2+H*KAPPA)/2.D0
  185.       A = DSQRT(A2)
  186.       B1 = H**2/8. -XA/2.
  187.       K2 = ( A +B1 -XA )/(2.*A)
  188.       KC2 = (A -B1 +XA )/(2.*A)
  189.       KC = DSQRT(KC2)
  190.       AP = A+XA-XP
  191.       AM = A-XA+XP
  192. C     IF (DABS(AM/A).LT..001) AM = XP+XA*(DSQRT(3.D0+H/K(3))-1.)
  193.       XAP = XA-XP
  194.       ALPHA = AM/AP
  195.       FP = AP**2/(4.*A*XAP)
  196.       RT = DSQRT(K2+AM**2/(4.*A*XAP))
  197.       RA = DSQRT(A)
  198.       IF (K2.GT.0) GO TO 19
  199.       AK = .0D0
  200.       KK = PI/2.D0
  201.       EK = KK
  202.       PIK = KK/DSQRT(FP)
  203.       GO TO 17
  204.  19   AK = DSQRT(K2)
  205.       CALL CEL1 ( KK, AK, IER )
  206.       CALL CEL2 ( EK, AK, 1.D0, KC2, IER )
  207.       CALL CEL3 ( PIK, KC, K2, FP )
  208.  17   FJ1 = -KK/RA
  209.       AL2M = -(4.*A*XAP)/AP**2
  210.       A2OM = -AM**2/(4.*A*XAP)
  211.       AKK = KC2*ALPHA**2 +K2
  212.       IF (DABS(ALPHA).LT..01) GO TO 25
  213.       FJ2 = (AP*PIK/(2.*XAP) -KK)/(RA*AM*RP)
  214.       FJ3 = FJ2/(AM*RP) -(A*PIK/(2.*XAP**2) -4.*A**2*(A2OM*EK
  215.      * -AKK*KK/AL2M +(A2OM**2-K2)*PIK)/(ALPHA*AP**3*AKK))
  216.      * /(RA*AM*RP**2)
  217.       GO TO 30
  218.  25   FJ2 = KK/(2.*RA*RP*XAP)
  219.       IF (K2.GE..0001) TEMP = (KK-EK)/K2
  220.       IF (K2.LT..0001) TEMP = PI/4.
  221.       FJ3 = (.5D0*TEMP-KK)/(2.*RA*(RP*XAP)**2)
  222.  30   S(1) = YP*FJ2 -H2*Y*FJ3
  223.       S(2) = 18.D0*FJ1/RP**2 -Y*FJ2
  224.       S(3) = (AP*KK/RA -2.*RA*EK)*RP**2/3.D0
  225.  5    IF (DABS(E-EMID).GT.DDE.OR.DLE.LT.0.) GO TO 20
  226. C
  227. C     CORRECT FOR DISCONTINUITY AT E = EMID
  228.       NUM(2) = -Q/KAPRM**2
  229.       NUM(1) = R1/KAPRM**2 +H2*NUM(2)/KAPRM
  230.       NUM(3) = RP/3.D0
  231.       G = 1.5707963D0*DABS(KAPRM)/(2.*K(3))
  232.       S11 = G*NUM(1)
  233.       S12 = G*NUM(2)
  234.       DO 8 L = 1, 3
  235.  8    S(L) = S(L) + G * NUM(L)
  236.       GO TO 20
  237. C
  238. C     CALCULATE OUTSIDE S = S11, S12
  239.  40   RT = DSQRT((KAPPA-H*XP)**2-4.*XP**3)
  240.       S12 = -PI*Y/(RP*RT)
  241.       IF (ISIG.GT.3) PRINT 5000, E, S12
  242.  5000 FORMAT (' E = ',E12.4,'  S12 = ',E13.6)
  243.       GO TO 20
  244. C
  245. C     CALCULATE COEFFICIENT OF LN TERM,E NEAR EL
  246.  50   XC = (H/3.D0)**2
  247.       P = KAPRM - RP*XC
  248.       TEMP = -(DABS(P)/DSQRT(3.*XC))
  249. C     TEMP = TEMP*((2.+3.*EPS/DE)*ALOG(DABS(H*P*(DE+EPS)/(81.*XC**2)))
  250. C    * +(2.-3.*EPS/DE)*ALOG(DABS(H*F*(DE-EPS)/(81.*XC**2))) -6.)
  251.       P2 = P**2
  252.       S(3) = TEMP*RP/3.D0
  253.       S(2) = -TEMP*(Q+RQ*XC)/P2
  254.       S(1) = TEMP*R1/P2 +H2*S(2)/P
  255.       E = E + EPS
  256. C
  257. C     COMMON TERMINATION
  258.  20   IF (ISIG.GT.3) PRINT 3000, E, ( S(L), L= 1,3 )
  259.  3000 FORMAT (' E = ',E12.4,' S S = ',3E14.8 )
  260. C     IF (IABS(ISIG).GT.1) PRINT 4000, FJ1, FJ2, FJ3
  261.  4000 FORMAT (' FJ S ',3E14.8 )
  262.       RETURN
  263.       END
  264. C
  265.       REAL*8 FUNCTION DARCOS(X)
  266.       REAL*8 X,PI
  267.       COMMON/PI/PI
  268.       IF (X**2-.0001) 1, 1, 2
  269.  1    DARCOS = PI/2. - X - X**3/6.
  270.       RETURN
  271.  2    DARCOS = DATAN( DSQRT(1.D0-X**2)/X)
  272.       IF (X.LT.0.) DARCOS = DARCOS + PI
  273.       RETURN
  274.       END
  275. C
  276. C  .............................................................................
  277. C      SUBROUTINE CEL1
  278. C
  279. C      PURPOSE
  280. C         CALOULATE COMPLETE ELLIPTIC INTEGRAL OF FIRST KIND
  281. C
  282. C      USAGE
  283. C         CALL CEL1(RES,AK,IER)
  284. C
  285. C      DESCRIPTION OF PARAMETERS
  286. C         REC   - RESULT VALUE
  287. C         AK    - MODULUS (INPUT)
  288. C         IER   - RESULTANT ERROR CODE WHERE
  289. C                 IER=0  NO ERROR
  290. C                 IER=1  AK NOT IN RENGE -1 TO +1
  291. C
  292. C      REMARKS
  293. C         THE RESULT IS SET TO 1,E75 IF ABS(AK) GE 1
  294. C         THE RESULT IS SET TO 1,0E38 IF ABS(AK) GE 1
  295. C         FOR MODULUS AK AND COMFLEMENTARY MODULUS CK,
  296. C         EQUATION AK*AK*CK*CK=1,0 IS USED
  297. C         AK MUST BE IN THE RANGE -1 TO +1
  298. C
  299. C      SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  300. C         NONE
  301. C         SPMAX
  302. C
  303. C      METHOD
  304. C         DEFINITION
  305. C         CEL1(AK)=INTEGRAL(1/SQRT((1+T*T)*(1+(CKKT)**2)), SUMMED
  306. C         OVER T FROM 0 TO INFINITY),
  307. C         EQUIVALENT ARE THE DEFINITIONS
  308. C         CEL1(AK)=INTEGRAL(1/(COS(T)SQRT(1+(CK*TAN(T))**2)),SUMMED
  309. C         OVER T FROM 0 TO PI/2),
  310. C         CEL1(AK)=INTEGRAL(1/SQRT(1-(AK*SIN(T))**2),SUMMED OVER I
  311. C         FROM 0 TO PI/2), WHERE K=SQRT(1,-CK*CK),
  312. C         EVALUATION
  313. C         LANDENS TRANSFORMATION IS USED FOR CALCULATION,
  314. C         REFERENCE
  315. C         R,BULIRSCH, 'NUMERICAL CALOULATION OF ELLIPTIC INTEGRALS
  316. C         AND ELLIPTIC FUNCTIONS',HANDBOOK SERIES SPECIAL FUNCTIONS,
  317. C         NUMERISCHE MATHEMATIK VOL, 7, 1965, PP, 78-90,
  318. C
  319. C   ............................................................................
  320. C
  321.       SUBROUTINE CEL1(RES,AK,IER)
  322.       IMPLICIT REAL*8 (A-H,O-Z)
  323.       IER=0
  324.       ARI=2.
  325.       GEO=(0.5D0-AK)+0.5D0
  326.       GEO=GEO+GEO*AK
  327.       RES=0.5D0
  328.       IF(GEO)1,2,4
  329.     1 IER=1
  330. C   2 RES=1.E75
  331.     2 RES = SPMAX(1)
  332.       RETURN
  333.     3 GEO=GEO*AARI
  334.     4 GEO=DSQRT(GEO)
  335.       GEO=GEO+GEO
  336.       AARI=ARI
  337.       ARI=ARI+GEO
  338.       RES=RES+RES
  339.       IF(GEO/AARI-0.9999)3,5,5
  340.     5 RES=RES/ARI*6.283185D0
  341.       RETURN
  342.       END
  343. C
  344. C     ..........................................................................
  345. C
  346. C        SUBROUTINE CEL2
  347. C
  348. C        PURPOSE
  349. C           COMPUTE THE GENERALIZED COMPLETE ELLIPTIC INTEGRAL OF
  350. C           SECOND KIND,
  351. C
  352. C        USAGE
  353. C           CALL CEL2(RES,AK,A,B,IER)
  354. C        DESCRIPTION OF PARAMETERS
  355. C           REC   - RESULT VALUE
  356. C           AK    - MODULUS (INRUT)
  357. C           A     - SONSTANT TERM IN NUMERATOR
  358. C           B     - FACTOR OF QUADRATIC TERM IN NUMERATOR
  359. C           IER   - RESULTANT ERROR CODE WHERE
  360. C                   IER=0  NO ERROR
  361. C                   IER=1  AK NOT IN RENGE -1 TO +1
  362. C
  363. C        REMARKS
  364. C           FOR ABS(AK) GE 1 THE RESULT IS SET TO 1,E75 IF B IS
  365. C           POSITIVE, TO -1,E75 IF B IS NEGATIVE,
  366. C           FOR ABS(AK) GE 1 THE RESULT IS SET TO 1,0E38 IF B IS
  367. C           POSITIVE, TO -1,0E38 IF B IS NEGATIVE,
  368. C           SPECIAL CASES ARE
  369. C           K(K) OBTAINED WITH A = 1, B = 1
  370. C           E(K) OBTAINED WITH A = 1, B = CK*CK WHERE CK IS
  371. C           COMPLEMENTARY MODULUS,
  372. C           B(K) OBTAINED WITH A = 1, B = 0
  373. C           D(K) OBTAINED WITH A = 0, B = 1
  374. C           WHERE K, E, B, D DEFINE SPECIAL CASES OF THE GENERALIZED
  375. C           COMPLETE ELLIPTIC INTEGRAL OF SEOOND KIND IN THE USUAL
  376. C           NOTATION, AND THE ARGUMENT K OF THESE FUNCTIONS MEANS
  377. C           THE MODULS,
  378. C
  379. C        SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED
  380. C           NONE
  381. C           SPMAX
  382. C
  383. C        METHOD
  384. C           DEFINITION
  385. C           RES=INTEGRAL((A+B*T*T)/(SQRT((1+T*T)*(1+(CK*T)**2))*(1+T*T))
  386. C        SUMMED OVER T FROM 0 TO INFINITY),
  387. C        EVALUATION
  388. C        LANDENS TRANSFORMATION IS USED FOR CALCULATION,
  389. C        REFERENCE
  390. C        R,BULIRCH, 'NUMERICAL CALCULATION OF ELLIPTIC INTEGRALS
  391. C        AND ELLIPTIC FUNCTIONS' , HANDBOOK SERIES SPECIAL FUNCTIONS,
  392. C        NUMERISOHE MATHEMATIK VOL, 7, 1965, PP, 78-90,
  393. C
  394. C     ..........................................................................
  395. C
  396.       SUBROUTINE CEL2(RES,AK,A,B,IER)
  397.       IMPLICIT REAL*8 (A-H,O-Z)
  398.       IER=0
  399.       ARI=2.
  400.       GEO=(0.5D0-AK)+0.5D0
  401.       GEO=GEO*AK + GEO
  402.       RES=A
  403.       A1=A+B
  404.       B0=B+B
  405.       IF(GEO)1,2,6
  406.     1 IER=1
  407.     2 IF(B)3,8,4
  408. C   3 RES=-1.E75
  409.     3 RES=-SPMAX(1)
  410.       RETURN
  411. C   4 RES=1.E75
  412.     4 RES=SPMAX(1)
  413.       RETURN
  414.     5 GEO=GEO*AARI
  415.     6 GEO=DSQRT(GEO)
  416.       GEO=GEO+GEO
  417.       AARI=ARI
  418.       ARI=ARI+GEO
  419.       B0=B0+RES*GEO
  420.       RES=A1
  421.       B0=B0+B0
  422.       A1=B0/ARI+A1
  423.       IF(GEO/AARI-0.9999)5,7,7
  424.     7 RES=A1/ARI
  425.       RES=RES+0.5707963D0*RES
  426.     8 RETURN
  427.       END
  428.       REAL*8 FUNCTION SPMAX(K)
  429.       SPMAX=1.0D38
  430.       RETURN
  431.       END
  432.